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FILTERING  INTERPOLATORS 


I.  INTRODUCTION 

In  many  electro-optical  imaging  systems,  the  most  serious  impediment  to  detecting  low-contrast 
moving  targets  is  spatially  varying  background  clutter.  Multiframe  clutter-suppression  processing  can 
often  provide  an  effective  solution  to  this  problem.  In  its  simplest  form,  this  processing  involves 
subtracting  a  second  image  from  a  first  in  a  frame-differencing  signal  processor  (FDSP),  in  order  to 
suppress  die  fixed  background  on  a  pixel-by-pixel  basis.  (The  emphasis  here  will  be  on  a  first-order 
FDSP,  but  die  interpolators  described  will  also  be  applicable  to  higher-order  FDSPs,  or  to  any  other 
signal  processor  that  compares  multiple  frames  of  data.)  Since  one  frame  must  almost  always  be 
resampled  before  it  is  compared  to  the  other,  interpolation  error  contributes  to  clutter  leakage  in  die 
difference  frame.  The  key  idea  of  this  report  is  to  apply  a  spectral  filter  to  the  unshifted  data  that  matches 
the  filtering  effect  of  applying  the  interpolator  to  the  data  that  are  shifted,  thereby  minimizing  the  effect 
of  interpolation  error  in  the  difference  frame.  Theoretical  analysis  indicates  that  the  result,  which  we  will 
call  a  filtering  interpolator,  gives  large  clutter  leakage  reductions  with  minimal  loss  of  target  signature. 
Tests  on  simulated  imagery  confirm  target-to-clutter  ratio  improvements  of  a  factor  of  five  or  more  for 
those  highly  structured  scenes  examined,  in  the  absence  of  other  sources  of  error.  (The  simulation 
assumes  a  space-based  IR  surveillance  sensor,  but  that  is  not  important  to  the  results  presented  here.) 

To  set  up  the  mathematical  problem,  let  x*  denote  the  value  at  u  =  n  of  data  taken  from  a  continuous 
function y(u),  (i.e.,  x„  =  fin)).  Consider  two  sets  of  equally  spaced  samples  {x„}  and  {x„+,},  as  shown  in 
Fig.  1,  the  second  set  being  taken  at  points  shifted  to  the  right  a  distance  s  from  the  first.  Since 
integral-sample  shifts  are  effected  simply  by  reindexing  the  data  set  and  do  not  require  interpolation,  we 
will  assume  0  s:  s  51  1 .  The  first  data  set  {x„}  will  be  called  the  shift  data,  because  it  must  be  processed 
to  find  an  interpolated  value;  die  second  set  {x,+,}  will  be  called  die  no-shift  data.  Interpolations  are 
always  done  between  die  two  central  points:  a  new  value  y,  must  be  generated  from  die  {xj  to  compare 
to  x,.  (The  set  {x„}  is  then  moved  over  one  sample  for  die  next  interpolation.)  The  new  value  y,  is  to  be 
given  by  a  local  formula  (as  opposed  to  a  global  spline  method  [1,2]): 

An 

y,  -  E  mjb.  (1) 

*— (An-i) 

where  the  Am(s),  the  interpolation  coefficients,  are  functions  of  die  shift. 

We  define  what  will  be  called  standard  interpolators  as  those  for  which  4,(0)  =  6(/i);  that  is, 
Ao(0)  =  1, 4,(0)  =  0  for  n  *  0.  When  the  shift  is  zero,  these  interpolators  simply  reproduce  die  data 
to  which  they  are  applied:  y,  =  x,  =  xa.  But,  for  arbitrary  data,  no  interpolator  will  be  perfect  for  s  #  0; 
there  will  be  some  error  when  y,  is  compared  to  x,.  Interpolators  tend  to  reproduce  some  frequency 
components  (usually  low)  of  data  with  high  fidelity,  but  do  less  well  on  other  frequency  components 
(usually  high).  Thus,  the  shifted  data  have  been  through  a  spectral  filter,  while  the  unshifted  data  have 
not.  The  resulting  spectral  mismatch  can  lead  to  significant  clutter  leakage  in  the  difference  frame.  While 
contemplating  this  problem,  one  of  us  [Stocker]  realized  that  if  the  right  spectral  filter  could  be  applied 
to  the  no-shift  data,  then  the  difference  between  y,  and  x,  would  be  reduced,  perhaps  dramatically.  Of 
course,  such  a  filter  would  also  reduce  high-frequency  detail. 
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Fig.  1  —  The  interpolation  problem:  the  solid  curve  repreaeoti  the  ontlog  signal  and 
the  shift  and  no-shift  data  sets  are  separated  by  the  shift  s.  An  interpolated  value  must 
be  found  from  the  shift  data  set  to  compare  to  the  no-shift  data  valuer,. 


In  this  report,  we  develop  parameterized  families  of  filtering  interpolators  that  permit  die  user  to 
choose  the  best  trade-off  between  clutter  suppression  and  high-frequency  retention  for  a  particular 
application.  One-dimensional  data  will  be  considered  explicitly;  for  two  dimensions,  die  interpolators  are 
applied  to  each  dimension  separately. 

For  our  purposes,  the  shift  s  is  assumed  to  be  known.  Usually,  s  must  be  determined  from  die  data 
itself,  a  problem  addressed  by  Kuglin  and  Hines  [3]  and  Schaum  and  McHugh  [4],  The  shift  can  also  be 
found  by  searching  for  the  value  of  s  that  minimizes  die  difference  frame  [5].  Our  ultimate  aim  is  to 
minimize  clutter  in  an  FDSP,  so  s  must  be  accurately  known.  A  small  error  fir  in  die  determination  of 
shift  will  lead  to  an  error  in  the  difference  frame  of  approximately  f'(u)6s. 

2.  A  PRESCRIPTION  FOR  INTERPOLATORS 

One  way  to  generate  the  A,  is  to  select  a  set  of  basis  functions  fjp)  that  are  thought  to  fit  the  data 
fairly  well  and  require  that  these  functions  be  interpolated  exactly.  That  is,  the  A*  are  selected  to  satisfy 
the  equations 


m 

E  /„<»kw  -/.to.  ® 

«— (AT2-1) 

with  m  =  1 , . . . ,  N.  Note  that  for  s  =  0,  Eq.  (2)  is  satisfied  by  AJi 0)  =  fi (n).  Solution  of  these  equations 
leads  to  standard  interpolators  for  any  desired  basis  functions. 

For  the  polynomial  basis  functions  fm(x)  =  the  resulting  A’s  give  standard  Lagrange 
interpolation  [1],  that  is,  the  data  are  fit  with  a  polynomial  of  order  N  -  1,  from  which  y,  is  calculated. 
The  coefficients  for  these  Lagrange  functions  for  N  data  points  (LF-N)  interpolators  are,  for  even  N, 


AM 


m 

n 

t—vn-n 

ixo 


(s-k) 

(n  -  *)- 


(3) 
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Another  possible  basis  for  interpolation  is  the  set  of  trigonometric  functions  used  by  foe  discrete 
Fourier  transform  (DFT).  For  N  =  4,  these  are 

m  - 
m  - 
m  « 
m « 

Substitution  of  these  functions  in  Eq.  (2)  yields  foe  DFT-4  interpolator  that  was  shown  by  Lucke  [6]  to 
give  better  performance  than  LF-4  when  foe  data  sampling  rate  is  low. 

Interpolators  derived  from  Eq.  (2)  are  of  foe  standard  type:  when  applied  with  s  =  0,  they  simply 
reproduce  foe  data.  We  now  rewrite  Eq.  (2)  to  include  foe  condition  that  foe  interpolator  be  applied  to 
foe  no-shift  data,  but  only  to  filter  it,  not  to  shift  it: 


1, 

cos(kx/2), 

sin(mc/2), 

cos(tx). 


m 


E 

m—i/n-n 


m 


E 

(W2-0 


fjfl  ♦  s)Am(0). 


(5) 


The  left-hand  side  of  Eq.  (5)  is,  as  before,  foe  result  of  applying  foe  interpolator  to  foe  shift  data. 
The  right-hand  side  is  foe  result  of  applying  it  to  foe  no-shift  data.  The  requirement  that  foe  two  be  equal 
preserves  foe  exact  interpolation  of  foe  basis  functions  in  foe  sense  that,  when  foe  difference  of  foe  two 
processed  data  sets  is  taken,  foe  result  is  zero.  (Filtering  interpolators,  when  applied  with  s  ■  0,  do  not 
reproduce  foe  data  to  which  they  are  applied;  that  is,  from  Eq.  (1),  y,  *  x,.)  Equation  (5)  allows  foe 
A„(s)  to  be  found  for  a  given  set  of  basis  functions  once  foe  4,(0)  are  chosen.  In  what  follows,  we  will 
give  forms  for  arbitrary  even  N  where  possible  but  will  usually  restrict  our  attention  to  N  =  4.  Results 
for  N  »  6  will  be  given  in  Section  7. 

We  will  require  that  foe  interpolator  reproduce  a  constant  signal  exactly  (i.e.,  /,(*)  =  1  so  that 
Ed.(j)  =  1  V  s),  and  that  foe  no-shift  filter  be  symmetric  about  n  =  0  so  that  its  phase  is  a  linear 
function  of  frequency.  The  4,(0)  are  thereby  determined  up  to  a  single  parameter  a: 


A.,(°)  *  a* 

4,(0)  -l-2o, 
4(0)  - 
4(0)  -  0. 


(6) 


Only  three  samples  of  foe  no-shift  data  set  have  nonzero  coefficients;  when  <2  =  0,  only  one  sample  has 
a  nonzero  coefficient  and  Eq.  (5)  reduces  to  Eq.  (2). 


When  foe  4(0)  from  Eq.  (6)  are  substituted  into  Eq.  (5),  along  with  foe  LF-4  or  DFT-4  basis 
functions,  foe  resulting  equations  can  be  solved  for  foe  Am(s).  The  interpolation  coefficients  for  LF-4  are 
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AJs)  =  I  [6a  -  2(1  ♦  3a)  s  +  3 s2  -  s*\, 
o 

Ms)  =  ^  [2(1  -  2a)  -  (1  -  6a)  s  -  2sa  ♦  s3], 

0) 

At(s)  =  i  [2a  +  2(1  -3a)  s  +s2  -  s*), 

A2(s)  =  1  [  -  (1  -  6a)  s  +  s*]. 

For  DFT-4  we,  find 

AJs)  -  4  [1  -  2(1  -  2a)&m(vs/2)  -  (1  -  4o)cos(«)], 

4 

A^s)  =  1  [1  +  2(1  -  2a)cos(«/2)  ♦  (1  -  4a)cos(«)J, 

4  (8) 

A,(s)  -  A  [1  ♦  2(1  -  2o)sin(xr/2)  -  (1  -  4o)cos(xs)], 

4 

Aj(j)  »  A  (1  -  2(1  -  2o)cos(xr/2)  ♦  (1  -  4o)cos(xr)]. 

4 

At  this  point,  we  introduce  the  interpolation  function  h(x)  defined  by  the  equation 

h(n  -  s)  -  A.(s),  (9) 

and  A(x)  =  0  for  |x|  >  Nf 2.  The  interpolation  function  hQc)  is  well-defined  because  0  <L  s  j£  1.  Note 
that  h(  1)  =  /*(— 1)  =  a  and  h( 0)  =1-2 a.  For  DFT-4,  the  interpolation  function  is 


h(x)  =  A  [1  ♦  2(1  -  2o)cos(xx/2)  ♦  (1  -  4a)cos(«)].  (10) 

4 

With  a  =  0,  this  is  the  standard  DFT-4  interpolation  function  given  by  Luclce  [6].  Figure  2  plots  die 
LF-4  and  DFT-4  interpolation  functions  for  a  =  0  (standard  interpolation)  and  for  a  =  0.21  and 
a  =  0.2,  respectively,  (die  best  values  for  filtering  interpolation,  as  we  see  in  the  next  section). 


3.  OPTIMIZATION 

The  value  of  die  parameter  “ a ”  remains  to  be  determined.  The  motivation  for  filtering  the  no-shift 
data  is  to  perform  the  same  spectral  filtering  that  die  interpolation  process  performs  on  the  shift  data, 
thereby  reducing  the  difference  between  the  two  data  sets.  In  principle,  this  means  matching  both  the 
amplitude  and  phase  of  the  digital  filters  applied  to  die  two  data  sets.  Here  we  will  concentrate  on 
matching  amplitudes,  because  this  leads  to  a  very  effective  solution  to  die  FDSP  clutter  suppression 
problem.  We  will  see  later  how  a  considerable  degree  of  phase  control  is  built  into  Eq.  (3)  and  show  in 
die  Appendix  how  the  problem  of  matching  phases  of  die  filters  for  the  two  data  sets  can  be  approached 
directly.  Detailed  examinations  of  amplitude  and  phase  matching  in  interpolators  and  in  differencing  filters 
are  developed  by  Schaum  [7]. 
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-2-10  1  2 


X 


-2-10  1  2 
X 

R|.  2  —  Interpolation  functions  for  LF-4  and  DFT-4,  in 
standard  form  (a  «  0)  and  beat-filtering  form  (a  =>  0.21 
and  a  *  0.2) 


We  see  from  Eq.  (5)  that  die  filtering  coefficients  applied  to  the  no-shift  data  are  the  interpolation 
coefficients  for  s  -  0,  while  those  applied  to  die  shift  data  can  be  for  any  value  of  s.  Therefore,  in  order 
to  make  both  sets  of  coefficients  perform  (nearly)  die  same  spectral  filtering,  we  seek  the  value  of  “a” 
in  the  interpolation  function  that  will  make  die  amplitude  of  spectral  filtering  as  independent  of  shift  as 
possible.  The  value  of  “a”  that  fulfills  this  condition  will  be  called  die  best-filtering  “ a 

Letting  where /is  the  frequency  in  cycles/sample,  signify  the  shift-dependent  spectral  filter 
for  which  the  tap  weights  are  die  A„  we  use  standard  methods  [8]  to  find  die  magnitude  of  H(f,s): 


W,s)\> 


m  ft- 1 

E  2£ 

•— (AT2-1) 


m-m 


E 

n— <W3-t) 


'W^to®ofi(2«n0. 


(11) 


Figures  3(a)  and  (b)  show  contour  plots  of  \Htf,s)\  for  LF-4  and  DFT-4,  respectively,  in  their  standard 
form  (a  -  0).  Both  can  be  seen  to  produce  substantial  filtering  of  high  frequencies  for  nonzero  shifts; 
more  to  the  point,  the  filtering  is  a  strong  function  of  shift. 

We  now  investigate  how  to  choose  “«”  so  that  \H(f,s)\  is  as  independent  of  s  as  possible.  This 
amounts  to  making  die  contours  of  |A/(f,j)|  resemble  straight  lines  parallel  to  the  shift  axis.  With  one 
parameter,  this  is  best  done  for  a  fixed  frequency  /  =  /;  since  most  imagery  data  are  dominated  by  low 
frequencies,  /  should  be  near  zero.  To  make  a  contour  line  flat,  we  will  require  that  die  derivatives  of 
the  function  with  respect  to  s  vanish  to  as  high  an  order  as  possible. 
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Fig.  3  —  Contour  plots  of  the  amplitude  of  H(f,s)  for  LF-4  and  DFT-4  with  a  -  0  and  a  equals  best-filtering  value.  For  the 
latter,  the  contours  are  remarkably  close  to  bring  straight  lines  parallel  to  the  shift  axis,  except  for  frequencies  near /  =  1/2. 


0.21,  and  for  DFT-4,  it  is  a  =  0.201  *  0.2,  which 


fori  =  1/2,  f,  small.  (12) 

In  Eq.  (12),  "/,  small”  for  DFT-4  means  that  cos(2inq£)  is  adequately  approximated  by  1  -  (2x7n£)2/2; 
then  the  solution  to  Eq.  (12)  is  independent  of  /.  For  LF-4,  the  expansion  of  die  cosine  must  be  carried 
to  fourth  order  in/  it  turns  (Hit  that  the  sum  of  all  die  second  order  terms  is  zero,  regardless  of  the  value 
of  “a.”  Vanishing  of  die  second  derivative  is  required  in  Eq.  (12)  because  symmetry  of  the  functional 
dependence  on  s  means  that  the  first  derivative  is  necessarily  zero  at  s  -  1/2.  Figures  3(c)  and  (d)  show 
that  the  contour  lines  of  \H(f,s)\  are  now  almost  (but  not  exactly)  straight  lines,  except  for  high 
frequencies.  That  is,  |//(f,s)|  is  almost  completely  independent  of  shift  for  the  low  frequencies  that 
dominate  most  real-world  image  data.  This  is  more  true  for  DFT-4  than  for  LF-4;  it  should  therefore  be 
a  better  interpolator  for  data  that  contain  more  high  frequencies,  or  equivalently,  for  data  that  are  more 
sparsely  sampled. 


For  LF-4,  the  desired  value  is  a  =  0.208  * 
were  found  by  requiring  that 

il  | -  o 
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If  the  contour  lines  are  very  close  to  straight  for  most  frequencies,  we  can  reasonably  expect  the  same 
parameter  value  to  work  for  a  wide  variety  of  applications,  because  the  functioning  of  the  interpolator 
will  be  nearly  independent  of  the  frequency  content  of  the  data  and  of  the  shift.  However,  for  data  that 
are  dominated  by  high  frequencies,  it  may  be  desirable  to  make  the  high-frequency  contours  straighter 
than  shown  in  Figs.  3(c)  and  (d).  This  can  be  done  by  increasing  “a”  slightly,  to  about  0.22  or  0.24, 
at  the  cost  of  having  contours  that  are  not  straight  (slightly  convex  upward)  at  lower  frequencies.  In 
Sections  6  and  7,  we  will  see  examples  of  interpolators  that  cannot,  with  any  choice  of  parameters, 
achieve  the  degree  of  contour  line  straightness  that  LF-4  and  DFT-4  do.  With  such  nonstraight  contours, 
parameters  may  need  to  be  adjusted  for  different  applications,  i.e.,  for  variable  frequency  content  of  the 
data  and  for  different  shifts. 

That  filtering  interpolators  can  be  constructed  with  nearly  straight  contours  is  of  particular  interest 
for  higher-order  frame  differencing.  A  second-order  frame  differencer  might,  for  example,  start  with 
three  frames,  subtract  the  second  from  the  first,  the  third  from  the  second,  and  then  difference  the  two 
resulting  first-differenced  frames.  If,  as  is  usually  the  case,  the  three  frames  are  not  coaligned, 
interpolations  must  be  performed  at  each  stage  of  the  process.  Clutter  leakage  can  be  minimized  by 
initially  arranging  the  frames  in  order  of  increasing  shift  (recall  that  shift  is  always  between  zero  and  one; 
this  may  require  reindexing  the  samples  of  one  or  more  frames),  applying  a  filtering  interpolator  with 
zero  shift  to  the  right-most  (cf,  Fig.  1),  and  applying  the  same  interpolator  with  appropriate  shifts  to  the 
other  two.  With  all  three  frames  thus  resampled  on  a  common  grid  and  nearly  identical  spectral  filtering 
applied  to  each,  the  second-order  difference  can  be  generated  without  further  concerns  about  registration 
and  with  minimal  effects  from  interpolation  error.  Obviously,  the  same  principle  applies  to  any  number 
of  frames  or  order  of  differencing,  as  well  as  to  any  other  signal  processor  that  compares  multiple  frames 
of  data. 

The  assumption  that  the  right-most  data  are  the  no-shift  data  is,  of  course,  purely  a  matter  of 
convenience.  The  left-most  data  could  equally  well  be  chosen:  the  sole  effect  of  this  choice  would  be  to 
reverse  the  order  of  the  tap  weights  applied  to  the  shift  data.  Thus,  if  several  sets  of  data  are  involved, 
any  one  can  be  taken  to  be  the  unshifted  data,  with  the  filtering  interpolators  derived  here  applied  directly 
to  those  shifted  to  the  left  and  applied  with  reversed  tap  weights  to  those  shifted  to  the  right. 

The  previous  discussion  has  concentrated  on  the  amplitude  of  H(f,s)  rather  than  its  phase.  The  effect 
of  the  phase  of  a  digital  filter  is  to  shift  the  Fourier  components  of  the  signal:  if  the  phase  of  H(fus)  is 
0.2  rad,  a  sinusoid  of  frequency /,  is  shifted  by  0.2/2*  cycles.  Thus,  the  phase  of  the  filter  and  the  shift 
of  the  signal  are  equivalent  concepts,  and  the  reason  that  we  need  not  consider  the  phase  of  H(f,s) 
explicitly  is  that  its  phase  requirements  are  built  into  Eq.  (S)  in  the  form  of  the  requirement  that  basis 
functions  be  exactly  shifted.  This  point  will  not  be  pursued  further  here,  except  to  exhibit  Figs.  4(a) 
through  (d),  which  show  contour  plots  of  the  phase  errors  for  LF-4  and  DFT-4.  The  phase  error  is  the 
difference  between  the  actual  phase  of  the  filter  and  the  phase  of  a  perfect  shifter  (e.g.,  the  Nyquist 
reconstruction  formula).  The  phase  of  a  perfect  shifter  is  a  linear  ftinction  of  frequency  with  slope 
2ts:  =  2 irsf.  Then  each  frequency  component  of  the  signal  is  shifted  by  a  distance  (<f>/2x  cycles)  x 
(l//per  cycle)  =  s,  i.e.,  shift  is  independent  of  frequency.  In  Fig.  4,  the  phase  error  is  always  zero  at 
a  shift  of  0.5,  because  of  the  symmetry  of  the  tapweights  [9],  which  is  /4_j(0.5)  =  ^2(0.5),  ^,(0.5)  = 
v*,(0.5).  In  Figs.  4(b)  and  (d),  the  phase  error  is  zero  at  a  frequency  of  0.25  because  DFT-4  interpolates 
this  frequency  exactly. 

4.  TESTING  THE  INTERPOLATORS 

As  a  first  step  to  seeing  how  well  the  interpolators  work  in  practice,  they  were  applied  to  a 
representative  point  spread  function  (PSF).  A  PSF  is  used  as  a  test  waveform  because  it  is  the  most 
rapidly  varying  signal  that  will  occur  in  most  data,  and,  therefore,  the  hardest  to  interpolate.  (For  an 
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SHIFT  SH|FT 

(c)  (d) 

Fig.  4  —  Contour  plots  of  the  difference,  in  degrees,  between  the  phsse  of  and  the  phase 
of  a  perfect  shifter,  for  LF-4  and  DFT-4 


example  of  data  that  vary  more  rapidly,  consider  a  scene  containing  a  positive-contrast  point  source 
located  very  close  to  a  negative-contrast  point  source;  fortunately,  such  occurrences  are  rare.)  The 
representative  PSF  chosen  is  a  unit-height  Gaussian,  with  standard  deviation  a  given  in  units  of  the 
sample  spacing.  By  changing  a,  we  can  see  how  well  an  interpolator  works  as  a  function  of  sampling 
rate. 


Figures  5(a)  and  (b)  show  results  in  which  “worst-phased”  sampling  (samples  miss  the  peak  of  the 
curve  by  the  maximum  amount)  is  assumed  because  it  is  bound  to  occur  at  many  points  in  real  data.  A 
standard  deviation  of  a  -  0.8  is  shown  because  it  is  a  stressing  case:  it  is  about  die  lowest  sampling  rate 
for  which  die  interpolators  perform  well  for  any  value  of  “a.”  Also  with  slighdy  less  than  two  samples 
per  foil  width,  half  maximum  (FWHM)  of  the  PSF,  it  is  about  as  low  a  sampling  rate  as  a  reasonable 
system  should  have  (a  *  1  or  higher  gives  substantially  smaller  interpolation  errors).  For  a  =  0,  the 
long-dashed  curve  of  Fig.  S  (the  result  of  applying  the  interpolator  to  the  no-shift  data)  is  coincident  with 
the  solid  curve,  and  the  difference  between  it  and  the  short-dashed  curve  (the  result  of  applying  the 
interpolator  to  the  shift  data)  is  substantial.  For  best-filtering  “a,”  the  maximum  difference  between  the 
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SAMPLE  NUMBER  SAMPLE  NUMBER 

(»)  (b) 

Fig.  S  —  Interpolated  curves  for  LF-4  (a)  and  DFT-4  (b).  The  solid  curves  represent  Gaussian  PSFs,  which  have  been  sampled 
at  the  +’».  The  short-dashed  curves  show  the  values  interpolated  from  these  samples.  The  iong-dasKed  curves  show  the  effect 
of  the  interpolator  on  the  no-shift  data,  for  all  possible  no-shift  data  sets  (i.e.,  on  the  solid  curve).  For  a  <■  0,  the  long-dashed 
curve  is  coincident  with  the  solid  curve. 

long-dashed  and  short-dashed  curves  is  reduced  by  about  an  order  of  magnitude;  a  =  0.1  is  an 
intermediate  case.  For  best-phased  sampling  (a  sample  occurs  at  die  peak  of  the  curve),  fits  are  much 
better  than  those  of  Fig.  5  for  a  *  0  but  not  noticeably  different  for  best-filtering  “a.”  For  a  1  and 
best-filtering  “ a ,”  the  two  dashed  curves  are  nearly  identical,  regardless  of  sample  phasing. 

Direct  testing  of  interpolators  on  real  data  is  very  desirable,  but  seldom  possible:  an  interpolated 
value  can  be  calculated,  but  die  true  value  to  which  it  should  be  compared  is  usually  not  available. 
However,  the  relative  performance  of  different  interpolators  can  be  compared  by  resampling  and 
differencing  simulated  imagery  and  calculating  the  root  mean  square  (r ms)  clutter  levels  of  the  difference 
frames.  With  computer-simulated  imagery,  we  have  die  luxury  of  studying  interpolation  error  in  isolation, 
separated  from  noise,  hardware  problems  such  as  line-of-sight  jitter,  or  other  complications.  Tests  on 
imagery  of  the  Earth’s  surface  [10]  have  verified  that  the  best-filtering  a’s  found  above  are  also  die  best 
values  to  use  on  imagery.  Tests  also  show  that  target-to-clutter  ratio  improvements  of  a  factor  of  five 
(14  dB)  or  more  are  obtainable,  compared  to  a  =  0  on  highly  structured  scenes,  including  an  areal  view 
of  a  city.  Improvements  on  low-structure  scenes  are  similar,  but  this  is  probably  not  important  because 
the  dominant  source  of  leakage  through  an  FDSP  for  these  scenes  is  noise,  not  interpolation  error. 

But  Figs.  3  and  3  show  that  the  price  of  this  additional  clutter  suppression  is  a  loss  of  high-frequency 
detail.  For  a  -  0,  only  the  shift  data  are  filtered  and  the  loss  of  detail  is  minimized.  For  best-filtering 
“ a ,”  both  data  sets  are  filtered,  and  the  loss  of  high  frequencies  can  be  substantial.  Equations  (7)  and 
(8)  constitute  one-parameter  families  of  interpolators:  the  user  can  choose  the  parameter  “a”  anywhere 
between  0  (minimum  high-frequency  loss,  maximum  clutter  leakage)  and  its  best-filtering  value  (greatest 
high-frequency  loss,  least  clutter  leakage),  depending  on  die  particular  target  strengths  and  clutter  and 
noise  levels  of  interest.  If,  for  example,  a  =  0. 1  reduces  clutter  leakage  below  die  level  of  sensor  noise, 
then  die  primary  effect  of  further  increasing  “a”  wi»i  be  to  reduce  target  signature. 
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5.  THEORETICAL  PERFORMANCE  ANALYSIS 

There  are  two  reasons  why  filtering  interpolators  work  better  than  standard  interpolators  in  preparing 
data  for  an  FDSP.  As  emphasized  above,  the  first  is  that  similar  spectral  filtering  is  done  on  both  data 
sets  (standard  interpolators  operate  on  only  one  data  set).  The  second  is  that  filtering  interpolators,  more 
than  standard  interpolators,  reduce  the  data’s  high  frequency  content.  This  includes  the  out-of-band 
contributions  from  the  displaced  replicas  of  the  spectrum  of  die  original  analog  signal  that  exist  in  the 
spectrum  of  the  sampled  data  and  are  aliased  into  the  passband  by  the  resampling  process.  This  latter 
point  is  addressed  in  the  Appendix,  which  gives  a  general  analysis  of  FDSP  performance  for  an  arbitrary 
interpolation  function,  or  kernel,  and  shows  how  the  kernel  acts  as  a  filter.  Figure  6  shows  the  frequency 
responses  of  the  standard  a*.1  best-filtering  interpolation  kernels  for  LF-4  and  DFT-4  (the  plots  were 
obtained  by  Fourier  transforming  the  relevant  four-point  kernels,  which  are  shown  in  Fig.  2). 

Figure  6  shows  that  the  reductions  in  the  peak  sidelobe  level  obtained  with  the  best-filtering  kernels 
are  large:  IS  dB  for  LF-4  and  30  dB  for  DFT-4.  Differences  in  stopband  behavior  between  the  two  types 
of  kernels  are  also  instructive.  The  LF-4  kernels  have  comparatively  high  sidelobes  with  wide  nulls 
around  integer  multiples  of  the  sampling  frequency  (i.e. ,  covering  the  repetitions  of  low  frequencies).  The 
DFT-4  kernels  have  much  narrower  nulls  but  also  have  a  much  lower  average  level  for  the  stopband 
sidelobes.  On  this  basis,  one  would  expect  LF-4  to  perform  well  on  low-frequency  data,  while  DFT-4 
would  be  a  better  choice  for  data  with  greater  high-frequency  content. 


Fig.  6  —  The  frequency  responses  of  LF-4  and  DFT-4  for 
a  -  0  and  a  equals  best-filtering  value.  The  plotted 
quantity  is  twenty  times  the  logarithm  of  the  magnitude  of 
the  Fourier  transform  of  the  interpolation  function. 
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Of  course,  the  cost  of  the  stopband  improvement  obtained  with  filtering  interpolators  is  additional 
in-band  attenuation  of  the  target  signal,  but  the  tradeoff  is  very  beneficial  in  terms  of  signal-to-clutter 
ratio  (SCR).  The  Appendix  shows  how  the  clutter  power  leakage  and  point  target  response  of  an  FDSP 
can  be  expressed  as  integrals  over  frequency  of  the  scene  and  point  response  power  spectra,  combined 
with  the  frequency  response  of  the  interpolation  process  (Eqs.  (A-8)  and  (A- 10)).  These  results  were  used 
to  calculate  interpolator  performances  for  the  specific  case  of  a  Gaussian  PSF  (a  =  0.8)  and  the 
representative  background  power  spectrum  given  in  Eq.  (A-9).  With  the  target  and  background  both 
scaled  to  unit  integrated  power,  the  difference  frame  SCR  given  in  Eq.  (A-ll)  is  a  good  measure  of  the 
signal  processing  gain  of  an  FDSP  for  different  interpolators,  compared  to  an  SCR  of  unity  for  the 
original  imagery.  Figure  7  shows  the  shift-dependent  SCR  gains  for  LF-4  and  DFT-4. 


Fig.  7  —  Relative  interpolation  performance  found  from  Eq.  (All)  of  LF-4  and  DFT-4 
with  a  *  0  and  a  equals  beat-filtering  value.  Note  that  wont  performance  ii  not 
neceaaarily  obtained  at  a  shift  of  0.5. 


The  gain  shown  in  Fig.  7  depends  on  two  factors:  the  clutter  reduction  in  die  difference  frame  and 
the  target  attenuation  due  to  no-shift  frame  filtering.  The  target  attenuation  penalties  incurred  with  the 
best-filtering  kernels  are  overwhelmed  by  the  dramatic  improvements  in  clutter  subtraction.  For  the  case 
considered,  SCR  improvements  of,  typically,  about  20  dB  (a  factor  of  10  in  amplitude)  are  obtained  by 
using  optimum  filtering  interpolators  in  lieu  of  standard  interpolators.  Note  that,  except  for  shifts  near 
zero  or  one,  the  best-filtering  DFT-4  kernel  yields  up  to  6  dB  more  gain  than  the  best-filtering  LF-4 
kernel  (for  the  model  spectrum  assumed).  Figure  7  also  shows  the  somewhat  counterintuitive  fact  that  a 
shift  of  1/2  is  not  always  the  worst  case  for  interpolation  error. 

The  reduction  in  target  signature  found  by  evaluating  the  integral  in  Eq.  (A-10)  with  and  without  the 
filtering  term  in  the  integrand  is  about  1  dB  for  best-filtering  “a,”  for  the  case  considered  here  (Gaussian 
PSF,  a  =  0.8).  This  result  is  for  integrated  signal  power.  For  many  systems,  the  reduction  in  peak  target 
signal  is  more  important;  it  can  be  seen  in  Fig.  5  to  be  about  2  dB,  or  a  factor  of  about  0.8  in  amplitude. 
The  maximum  possible  reduction  of  peak  target  signal,  which  occurs  for  a  target  that  has  only  one  sample 
on  it,  can  be  seen  from  Eq.  (6)  to  be  1  -  2a;  ordinarily,  data  are  sampled  often  enough  that  this  degree 
of  target  signature  reduction  will  not  occur. 
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The  typical  improvement  of  20  dB  given  above  is  better  than  the  14  dB  improvement  quoted  in 
Section  4  for  highly  structured  imagery,  partly  because  of  the  particular  model  background  power 
spectrum  used  here  and  partly  because  only  a  one-dimensional  (1-D)  problem  is  being  considered.  For 
isotropic  two-dimensional  (2-D)  data,  errors  resulting  from  interpolation  in  one  direction  are  statistically 
independent  of  errors  resulting  from  interpolation  in  the  other;  the  two  components  are  added  in 
quadrature  to  find  total  clutter.  Thus,  if  a  filtering  interpolator  reduces  one  component  of  clutter  by  some 
factor,  it  should  reduce  the  other  component,  and,  consequently,  the  total  clutter,  by  the  same  factor; 
i.e.,  clutter  reduction  is  independent  of  dimensionality.  But  target  signature  reduction  is  essentially 
multiplicative:  if  a  1-D  interpolation  reduces  a  target  signal  by  2  dB,  then  a  2-D  interpolation  will  reduce 
it  by  (about)  4  dB.  Thus,  for  properly  related  1-D  and  2-D  power  spectra  [11]  (i.e.,  both  are  obtained 
from  the  same  scene),  SCR  gain  will  be  lower  for  2-D  data  than  for  1-D  data. 

The  effect  of  the  interpolation  process  on  the  statistical  properties  of  the  noise  can  be  evaluated  by 
choosing  a  data  set  {*„}  that  consists  of  noise  with  zero  mean  and  variance  (x2)  =  a2.  It  is  easy  to  show 
from  Eq.  (1)  that  (y,)  =  0  and  that 


Stl  N- 1  Ntl-m 

lyj)  -  E  «H.!W  E 

m— (W2-1)  m-l  ■— <W2-1) 


(13) 


i.e.,  the  variance  of  the  interpolated  value  depends  on  the  autocorrelation  function  of  the  noise.  For 
uncorrelated  noise,  only  the  first  sum  contributes;  Fig.  8  plots  the  result. 


O 


Parameter  a 

Pig.  8  —  The  curvet  show  the  effect  of  the  interpolators  on  the  standard  deviation  of 
uncorrelated  noise,  at  a  function  of  the  parameter  "a,”  from  Eq.  (13).  For  5  =  0  (top 
curve),  LF-4  and  DFT-4  do  the  same  thing.  The  greatest  amount  of  noise  filtering, 
shown  in  the  lower  curves,  occurs  for  5  =  1/2. 
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6.  MORE  INTERPOLATORS/POSSIBLE  SPLINE  APPLICATION 

Equations  (7)  and  (8)  show  that  the  parts  of  the  An  that  depend  on  the  parameter  “a”  can  be 
separated  out.  Denoting  this  part,  which  is  the  difference  between  standard  ( a  =  0)  and  filtering  (a  *  0) 
interpolators,  by  AA„,  we  have 


AA_t  =  a  -  as, 

A 4,(s)  =  la  +  'has, 
AA,(s)  =  a  -  3as, 
AA2(s)  =  as. 


(14) 


for  LF-4.  These  AA’s  can  be  added  to  any  standard  N  =  4  interpolator  to  make  a  filtering  interpolator, 
and  the  same  can  be  said  for  the  A/t’s  obtained  from  Eq.  (8).  Interpolators  generated  in  this  manner  no 
longer  have  the  property  that  they  interpolate  a  simple  set  of  N-basis  functions  exactly  (which  is  also  true 
of  most  standard  interpolators).  Another  way  to  generate  new  interpolators  is  to  use  other  sets  of  basis 
functions  in  Eq.  5. 

For  example,  Eq.  (14)  can  be  added  to  the  coefficients  of  the  parameterized  cubic  interpolator 
described  by  Wolberg  [2]  to  yield 


A_t(s)  *  -4a  +  8a(l  +  s)  -  5a  (1  +  s)2  ♦  a  (1  *  s)3 

A0(s)  ■  1  -  (a  +  3)  s2  *  (a  ♦  2 )s3 

/i,(s)  =  1  -  (a  ♦  3)(1  -  s)2  «■  (a  +  2)(1  -  sf 

A2(s)  =  -4a  ♦  8a  (2  -  s)  -  5a  (2  -  s)2  *  a  (2  -  s)3 

We  take  the  parameter  a  (which  alters  the  frequency  response  of  the  interpolator)  to  be  - 1/2  (the  “image 
independent  optimal  choice”  of  Park  and  Schowengerdt  [12]).  For  this  interpolator,  it  turns  out  that  no 
choice  of  “a”  makes  the  contour  lines  of  |  H(f,s)  \ ,  shown  in  Fig.  9,  as  flat  as  those  for  LF-4  and  DFT-4, 
except  for  frequencies  very  near  /  =  0.  By  trial  and  error  on  simulated  imagery,  we  found  that 
(coincidentally)  a  =  0.2  works  best.  As  can  be  deduced  from  Fig.  9,  the  problem  with  this  interpolator 
is  that  it  works  well  (as  well  as  LF-4)  for  shifts  near  zero  and  one-half,  but  rather  poorly  (though  still 
better  than  its  standard  form)  for  shifts  near  one-quarter  or  three-quarters.  Other  values  of  a  and  “a” 
may  work  for  other  data,  but  the  user  must  proceed  at  his  own  risk. 

Another  candidate  for  modification  is  the  kernel  for  cubic  B-spline  interpolation  [2],  which  is  used 
here  as  a  local  rather  than  a  spline  interpolator.  (Pratt  [13]  also  uses  it  this  way  but  does  not  apply  it  to 
the  frame-differencing  problem.)  In  its  basic  form,  <4_,(0)  =  A,( 0)  =  1/6  and  A„(0)  =  2/3,  so  it  is  clearly 
a  filtering  interpolator;  in  fact,  it  is  the  first  one  we  used  to  test  the  idea.  The  interpolation  coefficients 
are 


*  a  -  as, 
-2a  +  3  as, 

♦  a  -  3 as, 

+  as. 


(15) 
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SHIFT 

Fig.  9  —  Contour  {riot  for  parameterized  cubic  (PC)  interpolation.  Moat 
contoun  have  the  same  value  near  *  *  1/2  as  near  5  =  0.  Hence  this 
interpolator  should  work  well  for  shifts  near  1/2,  but  not  for  shifts  near  1/4  or 
3/4. 


^(s)  *  I  [1  -  3s  ♦  3s2  -  s% 
o 

A,(s)  =  1  [4  -  6s2  +  3 s’), 

At(s)  -  i  [1  ♦  3s  ♦  3j*  -  3s3], 
6 

A2(s)  =  i  j\ 


(16) 


We  could  add  Eq.  (14)  to  these  to  change  their  filtering  properties,  but  inspection  of  Eqs.  (7)  shows  that 
Eqs.  (16)  are  equal  to  LF-4  with  a  =  1/6;  i.e.,  die  cubic  B-spline  interpolation  kernel  is  a  special  case 
of  parameterized  LF-4. 


The  DFT-4  interpolation  kernel  developed  here  can  also  be  applied  to  spline  interpolation.  A 
four-point  spline  interpolation  function  hQc)  should  be  (a)  everywhere  nonnegative  and  (b)  monotonically 
decreasing  to  zero  for  0  £  x  £  2  (cf.  Fig.  2  and  Ref.  [2]).  The  only  value  for  which  LF-4  has  these 
properties  is  a  =  1/6,  but  DFT-4  has  them  for  1/6  £  a  £  3/10.  We  suggest,  therefore,  that  DFT-4  with 
“ a ”  in  this  range  should  be  an  interesting  candidate  as  a  kernel  for  spline  interpolation. 


7.  HIGHER-ORDER  INTERPOLATORS 


For  some  applications,  higher-order  interpolators  may  be  desirable.  Lucke  [6]  gives  a  formula  for 
standard  DFT  interpolators  of  arbitrary  order,  and,  of  course,  die  standard  LF-N  interpolators  are 
available  in  many  software  packages.  Generating  filtering  interpolators  from  LF-N  or  DFT-N  for  N  >  4 
increases  rapidly  in  difficulty  because  of  the  number  of  equations  that  must  be  solved  and  die  number  of 
parameters  in  the  no-shift  filter.  Results  for  N  =  6  will  be  given. 
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For  six-point  interpolators,  imposing  the  same  conditions  as  for  Eqs.  (6),  the  4,(0)  are  given  by  two 
parameters,  a  and  b: 


AJO)  =  b,  4( 0)  -  a, 

^.,(0)  *  «,  4(°)  =  b, 

4,(0)  -  1  -  2a  -  2b,  4(0)  =  0. 


(17) 


With  two  parameters,  we  can  now  place  two  conditions  on  the  contours  of  \H(f,s) \ .  We  could,  for 
example,  extend  die  criterion  of  Eq.  (12)  to  the  fourth  derivative  to  obtain  maximal  straightness  of  the 
contour  lines  near  /  =  0  (condition  1).  Or  we  could  meet  the  criterion  of  Eq.  (12)  at  two  different 
frequencies  (condition  2). 


For  LF-6,  die  interpolation  coefficients  are 


-  2s*  +  s5}, 

♦  s4  -  t5}, 
(18) 

4(J)  =  m  -  5 (a  -  2 b))s  -  5[1  -  Ha  ♦  4b)]rs  ♦  s5}, 

A.2(s)  *  4(1  -  j),  4,(*)  =  4(1  -  s),  4(s)  =  4(1  -  s). 


4(s)  -  -1  [12a  «•  4[3  -  2(4<i  >  b)]s  *  8[1  -  3(o  ♦  4b)]r2  -  [7  -  20 (a  ♦  4b)]*3 
4(i)  =  _L  {24b  -  2[3  -4(4a  +  b)]s  -  {1  -  12(a  +  4b)]s2  ♦  [7  -  20 (a  ♦  4b)]*3 


With  this  interpolator,  as  for  the  parameterized  cubic  interpolator  of  Section  6,  the  parameters  that 
fulfill  condition  1,  namely  a  =  0.193,  b  =  0.101,  do  not  give  sufficiendy  straight  contours  for  larger 
/.  The  parameters  that  fulfill  condition  2,  with  f0  small  and/,  =  1/4,  namely  a  =  0.25S,  b  =  0.068,  do 
not  work  much  better  and  result  in  a  fairly  severe  loss  of  high  frequency  detail.  Once  again,  we  find  that 
no  choice  of  parameters  gives  really  straight  contours.  Hence,  for  this  interpolator,  we  should  expect  that 
different  applications  will  require  different  parameter  values.  Good  candidate  parameter  sets  to  start  with 
are  those  which  give  the  best  fits  for  our  Gaussian  PSFs,  a  =  0.24  and  b  =  0  (shown  in  Fig.  10(a))  and 
those  which  give  die  best  results  on  our  simulated  imagery,  a  =  0.26  and  b  =  0.05. 

For  DFT-6,  the  basis  functions  are 


/,(*)  - 1, 

/4(x)  -  cos(2xx/3). 

/2(x)  *  cos(xx/3), 

fM)  -  sin(2xx/3), 

(19) 

/,(x)  *  sin(xx/3), 

/#(x)  -  cos(xx). 

The  resulting  interpolation  function  is 


b(x)  -  4  11  ♦  2(1  -  a  -  3b)cos(xx/3)  ♦  2(1  -  3a  -  3b)cos(2«/3) 
6 

♦  (1  -  4a)cos(xx)] 


(20) 


15 


LUCRE  AND  STOCKER 


O) 


SAMPLE  NUMBER 


SAMPLE  NUMBER 


<») 


<b) 


Fig.  10  —  Interpolated  curve*  for  LF-6  (a)  and  DFT-6  (b). 
For  LF-6,  b  =  0  for  all  curve*. 


for  |jc|  ^  3  and  zero  otherwise,  and  A„(s)  =  h(n  -  s),n  =  -2,  -1,  0,  1,  2,  3.  The  parameter  values 
found  by  imposing  condition  1  are  a  =  0.236  *  0.24,  b  =  0.019  *  0.02.  In  this  case,  the  contour  lines 
of  | H(f,s) |  are  quite  straight  except  for  frequencies  close  to/=  1/2,  and  these  parameters  work  well  on 
die  Gaussian  PSFs,  as  shown  in  Fig.  10(b),  and  on  simulated  imagery.  For  a  =  0.8,  this  interpolator, 
in  its  best-filtering  form,  produces  errors  about  a  factor  of  two  smaller  than  DFT-4,  as  can  be  sera  by 
comparing  Fig.  10(b)  with  Fig.  5(b).  At  lower  sampling  rates,  DFT-6  maintains  a  performance  advantage 
over  DFT-4,  but  die  sampling  rate  cannot  decrease  much  before  die  data  become  so  poorly  sampled  that 
no  interpolator  works  well.  At  a  slightly  higher  sampling  rate  (a  ^  1),  die  performance  difference 
between  DFT-4  and  6  is  much  reduced:  for  well-sampled  data,  any  interpolator  gives  accurate  results. 

8.  CONCLUSIONS 

We  think  that  by  far  die  most  useful  interpolator  given  here  is  DFT-4  with  a  =  0.2  (best-filtering 
“a”).  We  recommend  it  for  general  applications,  even  if  die  resulting  interpolation  error  is  substantially 
smaller  than  the  noise.  LF-4  with  a  =  0.21  is  nearly  die  same,  but  produces  slighdy  more  high-frequency 
loss,  and,  as  can  be  sera  from  its  interpolation  function  (Fig.  2),  creates  (very)  small  discontinuities  in 
the  first  derivative  of  the  interpolated  data.  Only  if  minimizing  high  frequency  losses  is  a  major  concern 
is  it  likely  to  be  worthwhile  adjusting  “a”  to  a  smaller  value,  which  adjustment  (for  best  results)  must 
be  redone  for  data  sets  with  different  spectral  content.  Only  for  extremely  low-noise  data  would  a 
higher-order  interpolator  be  beneficial. 

If  the  best-filtering  value  cannot  be  used,  i.e.,  "a"  is  in  die  range  of,  say,  0  £  a  £  0. 15,  die  choice 
between  LF-4  and  DFT-4  is  determined  by  die  frequency  content  of  die  data:  die  difference  will  usually 
not  be  large,  but  LF-4  can  be  expected  to  give  better  results  on  data  that  are  more  strongly  dominated 
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by  low  frequencies  (or  are  more  densely  sampled).  Six-point  interpolators  will  improve  performance  in 
this  case,  but  we  have  scant  experience  on  which  to  base  recommendations.  Inspection  of  Fig.  10 
indicates  that  b  =  0,  a  £  0.1  should  work  fairly  well;  if  a  >  0.1  must  be  used,  probably  DFT-4  with 
a  =  0.2  would  be  a  better  (or  at  least  a  simpler)  choice. 

For  special  applications,  Section  6  gives  instructions  on  putting  the  interpolator  of  the  user’s  choice 
into  filtering  form.  The  filtering  interpolators  introduced  here  are  of  the  local,  convolutional  type  but  also 
lead  to  a  candidate  kernel  for  cubic  spline  interpolation  (Section  6  again).  The  new  interpolators  are  not 
appropriate  for  data  with  high  sampling  rates,  because  well-sampled  data  can  be  interpolated  with  high 
accuracy  by  standard  interpolators;  in  such  applications,  the  primary  effect  of  the  new  interpolators  would 
be  to  reduce  target  signatures.  But  for  moderate-to-low  sampling  rates,  we  have  found  SCR  improvements 
of  a  factor  of  five  or  more  for  point  targets  against  highly  structured  scenes,  when  interpolation  error  is 
the  dominant  source  of  clutter.  In  our  experience  with  simulations  of  a  state-of-the-art  IR  surveillance 
sensor,  clutter  leakage  due  to  interpolation  error  can  always  be  brought  below  the  level  resulting  from 
other  effects,  such  as  line-of-sight  jitter,  inaccuracies  in  determining  the  shift,  calibration  errors,  and 
detector  noise. 
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Appendix 

FOURIER  ANALYSIS  OF  INTRERPOLATORS 
FOR  FRAME  DIFFERENCING 


In  die  following  analysis,  a  continuous  signal  is  first  reconstructed  from  die  sampled  data  by 
convolving  with  an  interpolation  kernel;  the  result  is  then  resampled  at  die  specified  shift.  This 
interpretation  of  the  resampling  problem  will  show  explicitly  how  die  interpolation  function,  or  kernel, 
acts  as  a  filter  on  the  data.  Since  any  finite  interpolation  filter  has  a  frequency  response  of  infinite  extent, 
it  invariably  passes  out-of-band  components  in  die  Fourier  transform  of  die  sampled  data.  These 
out-of-band  components  alias  into  the  base-band  when  the  signal  is  resampled.  Expressions  for  the 
resulting  power  in  a  difference  frame  will  be  derived  for  a  representative  background  clutter  spectrum 
and  for  a  target  spectrum.  The  one-dimensional  case  is  done  here;  the  results  are  easily  extended  to  two 
dimensions. 

Let  y(x)  be  a  continuous  bandlimited  signal  with  Fourier  transform  Y{f)  defined  on  the  normalized 
frequency  band  -0.5  £  /  £  O.S.  This  signal  is  present  in  two  distinct  frames  and  is  sampled  at  equally 
spaced  but  differendy  phased  intervals  at  die  Nyquist  rate  (i.e.,  unity  sampling  interval),  as  shown  in 
Fig.  1.  The  discrete  frame  observations  (regarded  as  functions  that  are  defined  everywhere  but  are 
nonzero  only  at  the  sample  points)  and  their  Fourier  transforms  are  expressed  by 


/,(*)  -  yOO  E  «(*  -  m),  F,0)  *  £  Y(f  '  w)> 

m  m 

m  -  e  **  - « -  *).  m  -  E  r<f 

m  m 


(Al) 


where  m  ranges  from  -oo  to  o»,  and  x  is  the  shift.  Also, /,(x)  and  /2(x)  describe  die  shift  and  no-shift 
data  of  Fig.  1,  respectively. 

The  frames  are  to  be  resampled  to  spatially  align  the  signal  prior  to  differencing.  The  first  step  is 
linear  filtering  with  a  (continuous)  interpolation  kernel  with  finite  impulse  response  h(x)  (this  is  the 
interpolation  function,  cf.  Fig.  2)  and  frequency  response  H{f)  (the  magnitude  of  H(f)  for  LF-4  and 
DFT-4  is  shown  in  Fig.  6).  This  produces  a  pair  of  continuous  functions  denoted  by  g,  and  g2: 


g,(x)  -  /,(*)  *  A(x),  GW  *  H(J)  E  Y(f~  w>» 

at 

*,(*)  "/,(*)  *  h(x),  GJJ)  »  H(J)  £  Yif-  m**-. 

m 


(A2) 


To  obtain  a  pair  of  discrete  frames  that  are  spatially  co-aligned,  g,(x)  and  g2(x)  are  both  resampled 
at  x  *  n  +  s,  where  n  is  an  integer  (see  Fig.  1).  This  represents  a  shift  for  die  samples  from  which  gi(x) 
was  derived  but  not  for  die  samples  that  generated  g2(x).  The  resampled  discrete  frames  and  their  Fourier 
transforms  are  given  by 
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*„<*)  s),  G„(f)  =  £  J/(/-  n)  £  F(/-  m  - 

"  "  "  (A3) 

=  *2(*)  £  «<*  -  *  -  1),  <?»(/)  =  £«(/-»)  £  K(/-  rn  -  nye-fi**"*. 

ft  hot 

Ga  and  GM  are  defined  over  -  oo  <  /  <  oo ,  but,  since  they  are  die  Fourier  transforms  of  sampled  data, 
only  their  values  over  -0.5  ^  ^  0.5  are  of  interest.  Since  die  signal  spectrum  Y[f)  is  assumed  to  be 

bandlimited  to  this  range,  the  only  nonzero  terms  in  the  summations  over  m  are  m  =  -n,  when  -0.5 
£  /  £  0.5.  The  resampled  frame  transforms  then  simplify  to 


<W>  -*(/)£«(/-  -0.5  £  /  ^  0.5, 

<?*(/)  =  *(/)  £  «(/-»).  -0.5  i/  S!  0.5. 


(A4) 


The  spectral  content  of  the  difference  frame  is  therefore  given  by 

Daif)  =  GM(/)  -  G^CO  =  Y(f)  E(f,s),  (A5) 

where 

£(/,*)  =  £  /  -  «)[1  -  «**"]  (A6) 

ft 

is  a  shift-dependent  error  filter  defined  for  -0.5  £  /  £  0.5.  Thus  £(/»  is  a  weighted  sum  of  kernel 
response  samples  spaced  uniformly  at  intervals  of  the  sampling  frequency,  with  the  response  in  each 
repetition  interval  phase-shifted  by  an  amount  that  depends  on  the  frame  shift  s.  The  difference  frame  is 
always  zero  at  integer-valued  shifts.  Note  that,  for  die  ideal  sine  function  interpolator,  H(f)  =  Red <f) 
(i.e.,  H(f)  =  1  for  1/|  ^  1/2,  H(f)  -  0  otherwise),  so  that  E(f,s)  =  0  and  the  difference  frame  in 
Eq.  (A5)  is  identically  zero.  In  general,  however,  Eqs.  (A 5)  and  (A6)  show  that  the  difference  frame 
clutter  leakage  for  a  finite  interpolator  depends  on  the  out-of-band  frequency  response,  the  data  spectrum, 
and  the  shift. 


The  shift-dependent  power  due  to  clutter  in  the  difference  frame  is 


PJA  -  j”’  I *>„(/)!•  <//>  |n/)l!|£(/.t)l!  df.  (A7) 

This  expression  has  been  derived  for  an  arbitrary  deterministic  signal  with  Fourier  transform  Y(f).  For 
a  random  signal  (e.g.,  clutter  or  noise),  we  can  use  die  fad  that  die  expected  value  of  \Y(f)\7  over  die 
ensemble  is  the  discrete  signal  power  spectrum  S(j)  to  obtain  the  average  error  power  expression 


Pc  (*>  =  EPcto  =  12  Sif)\Eif,s)\7  df.  (A8) 

To  evaluate  Eq.  (A8),  we  need  a  representative  background  clutter  power  spectrum  to  use  for  S(f). 
We  have  found  the  following  formula,  which  gives  unit  integrated  power  over  the  range  -0.5  s;  /  £ 
0.5,  to  be  a  useful  model  for  imagery,  with  n  generally  in  die  range  of  4  to  6: 

S(f)  =  (n  ♦  1)(1  -  l\f\Y,  -0.5  £  /  s;  0.5.  <A9> 

Taking  n  =  4  yields  die  clutter  power  model  used  in  Eq.  (A8)  to  calculate  the  results  given  in  Section 
5. 
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The  target  response  attenuation  in  the  no-shift  frame  caused  by  the  use  of  a  filtering  interpolator  must 
also  be  calculated  to  obtain  the  signal -to-clutter  ratio  at  the  frame  differencing  signal  processor’s  output. 
Given  a  frequency-domain  description  T\f)  of  the  target  signal  shape,  the  output  target  power  in  the 
no-shift  frame  is  found  from  Eq.  (A4),  with  Y(j)  replaced  by  7\f),  to  be 

Pr  =  I!*,  |T(/)'2  IE  H(J  '  ^l2  df‘  (A10) 

J  -0.5  n 

The  output  target  power  in  the  shift  frame  could  be  similarly  calculated;  inspection  of  Eqs.  (A4)  and 
(A  10)  shows  that  it  is  shift-dependent.  We  now  assume  that  die  target  moves  sufficiendy  between  frames 
that  there  is  no  significant  interference  of  the  target  waveforms  when  the  difference  is  taken.  With  this 
assumption,  Eqs.  (A8)  and  (A10)  can  be  combined  to  compute  the  output  signal-to-clutter  ratio  given  in 
Fig.  7: 


SCR(s)  =  Ji-.  (All) 

Pc(s) 

The  results,  Eqs.  (A8),  (A10),  and  (All),  along  with  the  model  power  spectrum  Eq.  (A9),  were  used 
to  calculate  the  performance  figures  for  standard  and  filtering  interpolators  given  in  Section  S. 

We  have  attempted  to  improve  interpolation  performance  by  adjusting  the  parameter  “a:”  this 
determines  H(f)  and  the  error  fider  E(f,s).  Another  degree  of  freedom  can  be  developed  by  dropping  the 
condition  that  only  the  first  frame  is  shifted.  For  example,  we  might  “split  the  difference’’  between  the 
two  frames  by  moving  one  to  the  right  by  s/2,  the  other  to  the  left  by  s/2.  The  latter  is  the  same  as  a 
rightward  shift  by  1  -  s/2.  Inspection  of  Fig.  3  shows  that  |tf(/;s)|  is  the  same  function  of  frequency 
for  a  shift  of  1  -  s/2  as  it  is  for  s/2,  so  a  “two-way”  shift  produces  identical  (in  amplitude)  spectral 
filtering  of  the  two  data  sets.  But  it  does  not  necessarily  give  better  interpolation  results  than  a  one-way 
shift.  To  see  why  this  is  so,  we  first  note  that  if  resampling  is  not  to  be  done  on  one  of  the  original 
sampling  grids,  there  is  no  reason  why  it  must  be  done  at  s/2;  it  could  be  done  anywhere  in  the  [0,s] 
interval.  Such  a  two-way  shift  corresponds  to  replacing  x  -  n  -  s  as  the  argument  of  the  6  functions  on 
the  left-hand  side  of  Eq.  (A3)  by  jc  -  n  -  (3s,  with  0  £  /9  £  1. 


We  now  have  another  parameter  to  adjust  and  a  more  general  error  filter: 


£2(/,s)  =  Y  H(f  - 


(A12) 


=  2 j  Y  sin(xsn)[//(n  -  -  H(n  *  f)e'M1 ~V)M\, 

where  the  fact  that  H(J)  is  an  even  function  (it  is  the  Fourier  transform  of  h{x),  which  is  real  and  even) 
has  been  used.  Next,  we  rewrite  Eq.  (A6)  in  the  form 


£,(/*)  =  E  «(/-»)[!  -«*■"] 

m 

=  2 )  Y  sin (irsn)[H(n  -  -  H(n  *  f)e>m]. 

n*I 


(A  13) 
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The  relation  between  Eqs.  (A12)  and  (A  13)  is  hard  to  interpret  in  general,  but  die  special  case  0  = 
1/2,  s  =  1/2  is  instructive.  Equations  (A13)  and  (A12)  become,  respectively, 

£t(/,l/2)  *  2[H(l  -f)  *H(l  ♦/)],  (A14) 

and 

E2(f,  1/2)  -  2j[H(\  1  *f) ),  (A15) 

where  the  approximation  comes  from  die  fact  that  the  n  =  1  term  dominates  the  sum  (because  only  odd 
terms  give  nonzero  contributions,  which  decrease  rapidly  as  n  increases,  as  shown  in  Fig.  6). 

Equations  (A14)  and  (A15)  show  that  the  relative  performance  of  one-way  (EJ  and  two-way  (EJ 
interpolations  depends  on  the  relative  signs  and  magnitudes  of  H  just  above  and  below  unit  frequency. 
For  example,  the  relevant  sidelobes  for  DFT-4  (see  Fig.  6)  have  an  opposite  sign  before  the  absolute 
value  is  taken,  hence  tend  to  caned  in  Eq.  (A14)  and  reinforce  in  Eq.  (A15)  (i.e.,  |//(1  —  f)  +  H(  1  + 
/)|  <  |//(1  -  f)  -  H(  1  +  f)\  forO  SI  /  £  0.25).  Recalling  that  imagery  tends  to  be  dominated  by  low 
frequencies,  we  therefore  expect  that  a  one-way  shift  will  give  a  smaller  error  than  a  two-way  shift,  at 
least  for  0,s  -  1/2.  Trials  on  our  Gaussian  point  spread  functions  and  simulated  imagery  and  evaluation 
of  Eq.  (A1 1)  bear  out  this  expectation. 

In  general,  0  could  be  chosen  to  minimize  Eqs.  (A7)  or  (A8)  for  a  particular  interpolator,  signal 
power  spectrum,  and  shift.  We  have  not  investigated  this  point  extensivdy,  but  have  found  improvements 
in  clutter  suppression  for  “a”  near  zero  to  be  much  less  than  can  be  obtained  by  using  best-filtering  “ a .” 
For  best-filtering  "a,"  improvements  in  clutter  suppression  are  likdy  to  be  significant  only  for  very 
low-noise  data  (otherwise,  using  the  best-filtering  value  of  “a”  will,  by  itsdf,  reduce  interpolation  error 
bdow  the  noise  level).  For  DFT-4  with  best-filtering  “a,”  we  have  found  0  =  1  (conventional  one-way 
interpolation)  to  give  the  smallest  interpolation  errors.  For  LF-4  with  best-filtering  “a,”  a  two-way  shift 
with  0=1/2  tends  to  give  better  results  than  one-way  when  s  £  1/4. 
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